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ABSTRACT 


Increased  use  of  laminated  composites  has  pointed  out  the  need  for 
better  analytic  tools.  These  tools  must  be  able  to  correctly  account  for 

normal  shear  in  the  laminates  and  should  be  able  to  solve  nonlinear  problems. 

The  finite  element  method  can  be  applied  to  analyze  laminated  composites,  and 

this  research  presents  a  new  and  unique  method  to  include  nonlinear  effects 
and  normal  shear  effects.  The  finite  element  is  formulated  from  basic  elas¬ 
ticity  equations.  The  most  unique  characterisitc  of  the  element  is  the  manner 
in  which  normal  shear  is  handled.  The  normal  to  a  reference  surface  is 
allowed  to  not  only  rotate  but  also  to  change  shape.  Not  only  is  displacement 
continuity  imposed  at  lamina  interfaces,  but  also  slope  continuity  is 
guaranteed.  After  a  complete  general  formulation  of  the  finite  element  is 
presented,  the  method  is  specialized  to  plates  so  that  the  convergence 
characteristics,  the  accuracy,  and  the  applicability  of  the  element  can  be 
studied.  The  finite  element  does  an  excellent  job  of  analyzing  laminated 
composites,  although  there  is  some  degradation  in  accuracy  as  the  plates 
become  thinner.  This  finite  element  gives  the  analyst  the  ability  to  do 
nonlinear  analysis  of  laminated  composites  under  a  variety  of  loading  and 
boundary  conditions. 

A 


1.  Introduction 


PROBLEM  STATEMENT 

The  general  objective  of  this  dissertation  is  to  formulate  a  numerical 
procedure  for  the  large-displacement  analysis  of  laminated  composite  structures 
with  general  loading  and  boundary  conditions.  The  specific  numerical  procedure 
is  the  finite  element  technique.  The  translational  displacements  may  be  large, 
but  the  rotations  of  the  reference  surface  normal  will  be  infinitesimal. 

The  strains  will  remain  in  the  elastic  range.  To  examine  a  wider  range  of 
boundary  and  loading  conditions,  only  laminated  plates  will  be  considered. 
Several  loading  and  boundary  conditions  will  be  considered,  but  primary 
emphasis  will  be  on  normal  shear  behavior  and  those  conditions  which  have  a 
large  influence  on  shear. 

Laminated,  fibrous  composites  have  some  unique  advantages  over  isotropic 
materials,  but  they  also  present  some  additional  analysis  considerations. 

The  most  obvious  advantage  is  in  the  strength-to-weight  ratio  of  hi gh  strength 
composites  (1).  The  advantages  to  be  gained  through  weight  savings  are 
obvious.  The  other  advantage,  which  is  being  exploited  in  the  design  of  the 
forward  swept  wing  (2),  is  the  coupling  which  occurs  between  bending  and 
extension  in  unsymmetrically,  laminated  composites.  The  coupling  makes  it 
possible  to  adjust  the  stiffness  properties  of  the  material  to  meet  design 
limits  rather  than  only  adjusting  geometrical  sizes  of  the  structure. 

To  fully  capitalize  on  the  weight  advantage  of  composites,  the  designer 
should  optimize  the  size  of  the  structure.  Optimized  structures  are  sensitive 
to  collapse;  therefore,  a  complete  analysis  tool  should  include  at  least  geo¬ 
metrically  non-linear  analysis. 


A  unique  characteristic  of  composites  which  presents  a  challenge  to  the 
analyst  is  the  heterogeneous  character  of  the  fiber  and  matrix  constituents. 
Smearing  the  properties  of  these  two  constituents  and  treating  a  lamina  as 
an  orthotropic  material  in  the  plane  of  the  lamina  adequately  describe  the 
planar  response  of  the  lamina,  but  other  methods  must  be  used  to  handle  the 
out-of-plane  heterogeneity  of  laminates. 

These  other  methods  should  account  for  normal  shear  effects  since  these 
out-of-plane  effects  are  important  not  only  for  thick  structures,  but  also 
for  what  would  classically  be  considered  thin  structures  (3). 

An  analytical  method  must  include  nonlinear  geometrical  effects  and  normal 
shear  effects  if  laminated  composites  are  to  be  used  effectively.  By  including 
these  two  effects  in  a  generalized  computer  program,  structural  designers  and 
analysts  can  exploit  the  full  potential  of  composites. 

GENERAL  PROCEDURE 

The  procedure  to  be  followed  in  achieving  the  objective  of  this  research 
can  be  divided  into  four  steps.  First,  previous  work  in  composite  materials, 
nonlinear  geometrical  effects,  normal  shear  effects,  and  finite  element  analy¬ 
sis  is  examined.  This  area  is  covered  in  the  background  chapter.  Next,  the 
theory  and  methods  to  be  used  as  the  basis  for  a  numerical  analysis  of  the 
large  displacements  of  laminated  composite  structures  with  general  loading 
and  boundary  conditions  are  developed.  This  development  is  in  the  Theory 
chapter.  To  test  the  theory,  it  is  necessary  to  apply  it  to  specific  cases. 

For  this  dissertation  the  behavior  of  plates  under  a  variety  of  boundary  and 
loading  conditions  is  examined.  The  Results  chapter  presents  this  examination. 
Last,  the  degree  to  which  the  dissertation  objective  is  met  is  covered  in 
the  Conclusions  chapter. 


In  addition,  thfeei  appendices  are-  .Infiiudedt  The  first  appendix  presents 
the  derivation  of  the  strain-displacement  equations  for  the  unique  displace¬ 
ment  function  which  is  presented  in  this  work.  The  results  are  similar  to 
those  obtained  for  other  large-displacement  analyses,  but  are  derived  in  a 
unique  manner.  The  second  appendix  contains  the  derivations  of  the  stiffness 
matrix  for  a  constant  strain  triangular  plate  element  which  includes  not  only 
in  plane  translations  but  also  compatible  variations  in  displacements  through 
the  thickness.  This  material  is  related  to  the  information  in  the  Theory 
Chapter,  but  is  not  included  there  due  to  the  magnitude  of  the  manipulation 
required.  The  last  appendix  contains  an  explantion  and  derivation  of  the  cur¬ 
rent  stiffness  parameter  Vhi ch  is  used  in  the  nonlinear  analysis. 


In  this  chapter,  previous  work  in  composite  materials,  nonlinear  geo¬ 
metrical  effects,  normal  shear  effects,  and  finite  element  analysis  is 
examined  as  it  relates  to  this  research.  To  effectively  and  clearly  cover 
these  areas,  the  following  approach  is  used.  First,  the  methods  used  to  include 
three-dimensional  effects  in  plate  and  shell  theory  are  considered.  This 
consideration  will  concentrate  on  how  normal  shear  is  included  in  the  analysis 
of  isotropic  structures  and  1  ami'’ -ted  structures.  Next  the  methods  used  to 
include  these  three-dimensional  effects  in  finite  element  analysis  are  examined 
Last,  the  previous  work  in  non-linear  structural  analysis  is  reviewed. 

THREE-DIMENSIONAL  EFFECTS 

To  fully  appreciate  the  implications  of  three-dimensional  effects,  the 
basic  theories  must  be  understood.  A  recent  review  of  the  basic  theory  of 
shells  can  be  found  in  Gould's  book  (4).  One  of  the  earliest  works  on  ortho¬ 
tropic  shells  was  by  Hildebrand,  Reissner,  and  Thomas  (5).  This  paper  developed 
thin  shell  theory  from  the  three-dimensional  small  displacement  elasticity 
equations,  and  then  examined  the  effect  of  various  assumptions  on  the  formu¬ 
lation.  Laminated  plate  and  shell  theory  was  examined  extensively  in  Amburt- 
sumyan's  books(6,7).  Most  of  the  development  used  the  theory  of  thin  shells 
and  the  Kirchoff-Love  assumption.  The  last  chapter  does  contain  a  plane  stress 
type  of  development  in  which  transverse  shear  effects  are  included.  His 
book  on  anisotropic  plates  includes  a  large  section  on  analyzing  shear  in 
a  specific  case.  Books  by  Calcote  (8),  and  Jones  (9)  contain  reviews  of 


current  laminated,  composite,  plate  and  shell  theory.  The  basic  development 
in  both  books  still  includes  thin  plate  and  shell  theory  assumptioils  and  the 
Kirchoff-Love  assumption;  although  Jones  does  mention  that  transverse  stress 
effects  can  have  a  larger  impact  in  composite  laminates  than  is  experienced 
with  shells  made  of  isotropic  materials. 

To  examine  the  effect  of  transverse  shear  in  composites,  it  is  first 
necessary  to  look  at  What  was  done  with  isotropic  materials  to  include  the 
transverse  stress  effects.  One  of  the  earliest  works  was  presented  by 
Reissner  (10).  In  his  analysis  of  the  bending  plates,  the  normal  was  allowed 
to  rotate,  but  not  deform.  To  account  for  the  nonsatisfaction  of  equilibrium, 
shear  correction  factors  were  used.  In  1949,  Reissner  took  a  similar  approach 
in  analyzing  sandwich-type  shells  (11).  Another  shortcoming  of  earlier 
theories  was  accounted  for  by  Mindlin  (12)  when  he  included  the  effect  of 
transverse  inertia  terms  on  the  bending  of  plates. 

Another  way  to  account  for  transverse  effects  is  through  the  use  of  an 
oriented  continuum.  There  are  three  ways  to  formulate  oriented  continuum 
theories.  Various  authors  have  applied  Cosserat  surface  theory  to  shells  to 
include  couple  stresses  and  transverse  stress  effects  (13-16).  In  this 
theory  all  action  of  the  shell  is  simulated  by  a  surface  to  which  directors 
are  associated.  These  directors  cause  particles  to  not  only  have  trans¬ 
lational  displacements,  but  also  rotational  displacements.  A  second  directed 
continuum  theory  is  that  of  micromorphic  continua(17,18) . In  this  theory,  the 
assumption  is  made  that  the  continuum  is  composed  of  microelements  which  have 
their  own  microdisplacements  about  the  macroelement  center  of  mass  which  in 
addition  has  macrodisplacements  relative  to  the  global  system.  An  approach 
that  incorporates  elements  of  both  Cosserat  surface  theory  and  micromorphic 


RCAROOUCEj^N  US  AE A  COPIER 


[ 

I 

I 


continuum  theory  is  that  proposed  by  Mindlin  (19).  He  divided  the  continuum 
into  characteristic  parts  such  as  crystals,  grains  or  lamina.  The  motion  or 
the  entire  surface  was  described  by  the  displacement  of  the  center  of  mass  of 
these  characteristic  parts  and  microdisplacements  which  could  be  equivalently 
described  by  rotation  of  directors  assigned  to  the  center  of  mass. 

Although  there  is  experimental  evidence  that  this  type  of  continuum 
behavior  exists  when  an  clastic  layer  is  imbedded  between  two  more  rigid  layers 
and  experiences  a  shear  type  force  (20),  there  are  practical  problems  in 
obtaining  the  constitutive  equations.  Some  of  the  consti tui ti ve  properties 
are  represented  by  fifth  and  sixth  order  tensors  and  cannot  be  directly  obtained 
from  experiments.  The  usual  practice  to  obtain  these  material  properties  is 
to  obtain  the  strain  energy  density  for  an  equivalent  three  dimensional  problem 
and  equate  this  strain  energy  with  that  obtained  using  a  directed  continuum 
theory  (16).  Rather  than  use  a  directed  continuum  theory,  one  could  use  an 
equivalent  three-dimensional  problem  where  deformations,  which  vary  with  the 
thickness,  are  described  by  deformation  of  directors  attached  to  the  reference 
surface.  This  will  be  the  approach  followed  in  this  research. 

One  of  the  earliest  attempts  to  include  the  shear  stress  effects  for 
anisotropic  laminated  plates  was  made  by  Amburtsumyan  (7).  His  approach  was 
complicated  and  could  only  be  applied  to  plates.  Vang,  Norris,  and  Stavsky 
also  developed  a  theory  for  heterogeneous  plates  (21)  in  which  they  assumed 
that  the  in-plane  displacements  varied  linearly  with  the  thickness  coordinate. 
They  had  to  use  correction  factors  to  satisfy  boundary  conditions  and  equilib¬ 
rium.  The  next  significant  work  was  by  Pagano  (3).  He  analyzed  the  cylindri¬ 
cal  bending  of  a  semi-infinite  plate  using  elasticity.  His  work  has  been  used 
by  many  others  to  guage  the  accuracy  of  their  solutions.  Whitney  recently 
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applied  the  work  of  both  Amburtsumyan  and  Yang  to  specific  symmetric  and 
antisymmetric  laminated  composite  plate  problems  (22-28). 

The  next  logical  step  after  assuming  a  linear  variation  in  the  tangental 

displacement  through  the  thickness  was  to  assume  the  displacements  varied 

linearly  within  each  layer  of  a  laminate.  Grot  used  this  approach  for  lamin- 
•• 

ated  composites  and  divided  the  laminate  into  fiber  and  matrix  layers  (29). 
Each  layer  was  treated  as  isotropic  and  the  entire  formulation  was  worked  out 
in  terms  of  the  matrix  and  fiber  properties  and  displacements.  Srinivas  let 
the  displacements  vary  linearly  within  an  orthotropic  layer  composed  of  both 
fiber  and  matrix  (30)  and  used  his  formulation  to  analyze  the  vibrations  of 
laminated  plates.  An  approach  similar  to  that  developed  by  Srinivas  was 
used  to  formulate  the  equations  for  the  nonlinear  analysis  of  multilayered 
shells  (31).  In  this  paper,  the  similarities  between  this  approach  and  the 
results  obtained  using  directed  continuum  theories  were  also  examined.  Most 
recently  Pagano  developed  a  theory  using  Reissner's  variational  principle 
and  assumed  stress  fields  (32).  As  a  result  of  his  formulation,  he  obtained 
13N  field  equations  and  7N  edge  conditions  (where  N  is  the  number  of  layers) 
which  required  simultaneous  solution. 

All  of  these  approaches  except  the  most  recent  presented  by  Pagano  intro¬ 
duced  discontinuities  in  the  displacement  derivatives  at  the  layer  interfaces 
As  a  result  shear  correction  factors  had  to  be  used  to  guarantee  satisfaction 
of  the  equilibrium  equations.  Although  Pagano  did  guarantee  satisfaction  of 
equilibrium,  his  approach  results  in  a  large  system  of  linear  equations  which 
must  be  solved  simultaneously.  The  solution  of  this  system  would  be  even 
more  difficult  if  the  equations  were  nonlinear,  as  they  are  when  collapse 
is  analyzed.  This  initial  analytical  development  has  now  moved  into  finite 
element  analysis. 


FINITE  ELEMENT  ANALYSIS 


The  basic  theory  behind  finite  element  analysis  is  well  documented  (33-35) 
and  will  not  be  covered  here.  Rather,  the  methods  used  to  do  three-dimensional 
finite  element  analysis  will  be  examined. 

The  earliest  attempts  at  including  shear  in  finite  elements  were  in  two 
categories.  These  areas  were  sandwich  shell  analysis  (36)  and  three  dimensional 
continuum  analysis  (37).  The  sandwich  shell  type  elements  usually  analyzed 
the  outer  coverings  as  membranes  and  the  core  material  as  if  it  would  only 
carry  transverse  shear.  The  solid  hexahedra  is  the  usual  element  chosen  for 
three-dimensional  continuum  analysis  (34).  The  eight  node  trilinear  element 
and  the  twenty-one  node  serendipity  type  of  element  are  the  predominant  choice 
among  analysts  (35). 

Another  application  of  shear  in  finite  element  analysis  is  the  Mindlin 
bending  elements  (35).  In  these  elements  the  nodal  degrees  of  freedom  are  the 
transverse  displacement,  w,  and  two  rotations  of  the  normal  to  the  reference 
surface. 

Although  satisfactory  results  are  obtained  with  these  types  of  elements  for 
thick  structures,  they  could  not  be  used  to  analyze  thin  structures  without 
modification  (38).  The  problems  with  thin  structures  came  about  because  of 
the  phenomena  called  "locking"  in  which  the  shear  energy  predominated  rather 
than  going  to  zero  as  theoretically  required.  The  locking  also  occurred 
in  elements  which  used  linear  or  quadratic  shape  functions,  but  disappeared 
when  higher  order  shape  functions  were  used.  There  are  various  ways  of  reducing 
the  effect  of  locking  (39). 

The  simplest  method  which  has  been  presented  for  improving  thin  structure 
analysis  is  reduced  integration  (40,  41).  Other  methods  adjust  the  geometrical 
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properties  of  the  element  (39).  Recently,  penalty  functions  have  been  used 
successfully  (42).  The  result  of  all  methods  is  to  reduce  the  magnitude  of 
shear  stiffness  and,  hence,  significantly  reduce  the  magnitude  of  the  shear 
energy. 

Since  shear  is  more  dominant  in  anisotropic  materials  than  isotropic 
materials  (3),  it  was  logical  to  conclude  that  the  finite  elements  which  are 
used  for  composites  must  include  the  effects  of  transverse  shear.  The  first 
finite  elements  for  composites  which  included  shear  were  based  on  the  theory 
formulated  by  Amburtsumyan  (7)  and  applied  by  Whitney  (22).  This  theory  is 
a  Hindi  in  type  of  theory.  In  these  elements  (43,  44),  the  nodal  degrees  of 
freedom  were  the  three  translational  displacements  u,  v,  w,  and  two  rotations, 
and  which  describe  the  rotation  of  the  normal.  The  normal  for  these 
elements  remained  straight,  but  the  classical  Ki rchoff-Love  assumption  for 
normals  was  relaxed.  This  type  of  analysis  resulted  in  a  constant  value  of 
shear  per  layer.  This  meant  that  neither  compatibility  nor  equilibrium  were 
satisfied  at  the  lamina  interfaces. 

An  extension  of  this  theory  was  to  include  the  transverse  shears  at  each 
node  as  two  additional  degrees  of  freedom  (35,  45).  This  approach  did  give 
somewhat  more  accurate  results,  but  still  was  not  able  to  model  the  additional 
deformation  of  the  normal  which  occurs  in  a  laminated  structure  (45). 

The  next  extension  was  by  Mau,  Tong,  and  Pian  (46).  They  allowed  the 
normal  to  have  a  different  rotation  within  each  layer.  This  approach  has  also 
been  used  by  others  for  a  variety  of  problems  (47-49).  In  this  approach,  the 
shear  is  only  approximated  in  an  average  sense  per  layer  and  displacement 
continuity,  but  not  strain  compatibility  is  satisfied  at  lamina  interfaces. 

Since  shear  is  more  dominant  in  composites,  these  analyses  did  not  suffer 


from  locking  as  much  as  would  occur  for  isotropic  analysis,  but  the  phenomena 
still  existed.  The  effects  of  reduced  integration  and  higher  order  shape 
functions  was  studied  by  Reddy  (41).  He  showed  that  reduced  integration  and 
higher  order  shape  functions  did  improve  solution  accuracy. 

In  all  of  these  methods,  both  for  isotropic  and  anisotropic  materials, 
in  which  the  normal  is  allowed  to  rotate,  but  remains  straight,  a  shear 
correction  factor  had  to  be  used.  Several  authors  have  derived  formulas  for 
this  correction  factor,  which  adjusts  for  the  non-satisfaction  of  equilibrium 
by  the  constant  shears.  These  formulas  are  not  universally  applicable  (50), 
and  therefore,  the  confidence  is  low  for  the  usage  of  these  elements  in  varied 
loading  and  boundary  conditions. 

Another  recent  approach  for  laminated  finite  element  analysis  is  the 
hybrid  stress  element  (51).  This  element  uses  the  stresses  as  degrees  of 
freedom  along  with  the  usual  degrees  of  freedom.  It  also  employs  the  use  of 
Lagrangian  multipliers  to  ensure  the  satisfaction  of  equilibrium  at  layer 
interfaces.  It  then  uses  20  stress  degrees  of  freedom  per  layer  plus  dis¬ 
placement  degrees  of  freedom.  As  a  consequence,  the  number  of  required  compu¬ 
tations  is  extensive.  Although  one  stress  type  element  has  been  applied  to 
nonlinear  analysis  (52),  the  application  is  extremely  complicated  and  requires 
additional  computational  efforts. 

NONLINEAR  ANALYSIS 

As  Koiter  mentioned  in  a  paper,  linear  theory  cannot  adequately  describe 
the  collapse  behavior  of  structures  (53).  Since  solutions  of  nonlinear  equa¬ 
tions  are  not  readily  available,  numerical  methods  must  be  used.  In  this 
research,  the  finite  element  method  will  be  incorporated.  Good  reviews  of 
the  finite  element  method  and  the  application  to  nonlinear  shell  problems 


can  be  found  in  two  texts,  one  by  Brebbia  and  Connor  (54)  and  one  by 
Zienkiewicz  (33).  Brebbia  and  Connor  have  a  section  on  shell  analysis  which 
covers  the  inclusion  of  curvilinear  conditions  in  the  shape  function  and  the 
causes  of  the  problems  associated  with  rigid  body  displacements.  Zienkiewicz 
covers  the  solution  procedures  which  can  be  used  to  solve  the  nonlinear  shell 
equations.  Zienkiewicz  recently  collaborated  on  the  use  of  nonlinear  finite 
elements  to  study  one  dimensional  type  problem  (55).  A  larger  and  more 
complete  work  on  the  use  and  formulation  of  finite  elements  for  solving  non¬ 
linear  equations  was  recently  published  by  Dodds  (56).  Following  is  a 
synopsis  of  previous  research  in  nonlinear  finite  element  analysis. 

There  are  two  main  areas  to  consider  in  geometrically  nonlinear  structural 
analysis.  First,  is  the  formulation  of  the  nonlinear  equations.  The  second 
area  is  the  collection  of  methods  that  are  used  to  solve  these  equations. 

In  addition  to  these  two  main  areas,  the  work  that  has  applied  these  formula¬ 
tions  and  solution  technique  must  be  examined. 

All  formulations  of  nonlinear  structural  analysis  have  their  basis  in 
three-dimensional  continuum  mechanics  (57).  These  formulations  are  either 
built  with  a  Lagrangian  or  Eulerian  frame  of  reference.  In  the  Eulerian 
formulation,  the  equilibrium  equations  are  more  simply  expressed  (57).  In 
addition,  using  the  Eulerian  formulation,  nonlinear  problems  can  be  solved  as 
if  they  are  piecewise  linear.  Since  the  Lagrangian  formulation  is  based  on 
the  undeformed  configuration,  the  material  properties  in  the  undeformed  con¬ 
figuration  can  be  used  independent  of  the  deformations  (55).  This  property 
is  extremely  convenient  for  anisotropic  analysis.  There  are  several  deriva¬ 
tions  of  suitable  equations  for  shells  and  plates  based  on  the  three  dimen¬ 
sional  equations,  using  both  Lagrangian  and  Eulerian  systems;  these  are  also 


well  known,  and  detailed  derivations  can  be  found  in  elasticity  textbooks 
(58-60). 

There  have  been  some  applications  of  an  Eulerian  type  approach  for  solving 
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nonlinear  problems  using  finite  elements  (61).  The  disadvantage  of  this 
approach  is  that  coordinates,  strains,  stresses  and  material  properties  must 
be  transformed  at  each  load  step.  This  greatly  increases  the  number  of  required 
calculations,  especially  for  large  structures. 

The  Lagrangian  method  seems  to  have  become  the  most  favored  for  nonlinear 
finite  elements.  This  is  especially  true  for  anisotropic  and  laminated 
structures. 

The  various  methods  which  are  used  to  solve  nonlinear  equations  have  been 
studied  extensively  (57,  62-65).  These  include  the  following  methods:  incre¬ 
mental  stiffness,  predictor  corrector,  Newton-Raphson  iteration,  first-order 
and  second  order,  self-correcting  and  perturbation  (66).  Although  all  have 
been  used  successfully,  it  has  been  demonstrated  that  the  most  powerful  method 
are  the  Newton-Raphson  iteration  and  the  first-order  self-correcting  methods 
(66).  The  Newton-Raphson  method  usually  has  been  employed  due  to  its  simplicity. 

In  all  methods,  force  incrementation  has  been  used  in  the  pre-collapse 
portion  of  the  analysis  (57,65).  Post-collapse  analysis  has  usually  required 
displacement  incrementation  (66).  In  most  cases  collapse  has  been  determined 
by  examining  the  rank  of  the  structural  stiffness  matrix  (57,  64).  When  this 
matrix  becomes  singular  a  collapse  point  is  reached.  Recently,  Bergan  and 
others  proposed  the  use  of  a  current  stiffness  parameter  to  determine  collapse 
points  (64).  This  method  has  recently  been  successfully  applied  by  Noor  (65). 

The  current  stiffness  parameter  is  a  global  measure  of  the  structural  stiffness 
and  eliminates  the  need  for  determining  the  rank  of  the  large  finite  element 
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stiffness  matrix.  This  method  and  its  application  will  be  more  thoroughly 
explained  in  the  theory  chapter  and  Appendix  3. 

Nonlinear  analysis  has  been  applied  to  isotropic  materials  for  many  years. 
The  two  most  common  methods  have  been  determination  of  the  complete  nonlinear 
behavior  (67)  and  an  eigenvalue  analysis  to  determine  buckling  loads.  Nonlinear 
analysis  has  been  reviewed  in  many  articles  (57,  68,  69).  A  general  review  of 
buckling  and  collapse  analysis  can  be  found  in  the  text  by  Brush  and  Almroth 
(70).  Collapse  analysis  of  shells  is  covered  in  Gould's  text  (4).  Gould 
emphasizes  the  importance  of  the  effect  of  transverse  stresses  on  shell 
stabi 1 i ty . 

The  collapse  and  nonlinear  analysis  of  laminated  composite  structures  has 
only  been  studied  thoroughly  in  recent  years.  Calcote  discussed  the  stability 
analysis  of  laminated  composite  shells  (8),  but  his  analysis  used  the  Kirchoff- 
Love  assumption.  Khot  wrote  a  series  of  reports  which  analyzed  the  influence 
of  fiber  orientation,  non-homogenieties,  and  initial  imperfections  on  the 
buckling  of  laminated  composite  shells  (71-74).  More  recently,  Bauld  analyzed 
the  collapse  of  laminated,  composite  panels  using  finite  differences  (75). 
Missing  from  most  of  these  papers  is  experimental  data  to  compare  with  the 
analyses.  There  has  been  some  recently  published  experimental  data  (76). 

The  literature  on  complete  nonlinear  analysis  of  laminated  composites  is 
relatively  small  (48,  55,  77-81).  The  results  presented  have  been  of  complete 
plates  or  of  one  dimensional  type  structures.  Epstein  and  Glockner  have 
recently  presented  a  theory  which  includes  thickness  variations  in  transverse 
shear,  but  they  did  not  present  results  (80).  The  other  papers  used  either  a 
Mindlin  type  theory  (78,  81),  or  nonlinear  classical  laminated  plate  theory 
(77,  79). 
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SUMMARY 


Although  there  has  been  extensive  previous  work  in  composite  materials, 
nonlinear  geometrical  effects,  normal  shear  effects,  and  finite  element  analy 
sis,  there  has  not  been  a  great  deal  of  research  into  combining  these  four 
areas.  In  those  instances  where  these  areas  are  combined,  compatibility 
was  not  completely  satisfied  in  the  transverse  direction.  Presented  in  the 
following  section  is  the  theory  for  a  nonlinear  finite  element  for  the  analy¬ 
sis  of  laminated  composites  which  is  based  on  general  three  dimensional  con¬ 
tinuum  theory  and  does  completely  satisfy  compatibility  in  the  transverse 


direction. 


3.  Theory 


The  theory  behind  the  finite  element  solution  of  nonlinear  problems 
involving  laminated  fibrous  composites  can  best  be  understood  by  separating 
it  into  three  areas.  First  is  the  development  of  che  basic  equations  of 
strain  and  displacement  and  stress  and  strain.  Next  is  the  finite  element 
formulation.  Last  is  the  solution  procedure. 

BASIC  EQUATIONS 

There  are  several  criteria  which  must  be  considered  when  establishing 
relations  for  this  problem.  The  choice  of  these  criteria  is  influenced  by 
the  solution  method  .finite  elements,  the  type  of  problem,  geometrically 
nonlinear,  and  the  type  of  material,  laminated,  fibrous  composites.  Since 
a  geometrically  nonlinear  problem  is  studied,  the  strain-displacement  relations 
must  include  large  displacements.  As  has  been  shown,  accurate  analysis  of  com¬ 
posites  must  include  the  effects  of  normal  shear  (31).  To  include  these 
stresses,  the  displacements  must  be  dependent  on  the  thickness  coordinate. 

An  additional  consideration  is  the  reduction  of  the  adverse  effect  of  rigid 
body  displacements  in  satisfying  equilibrium  (82).  To  minimize  this  effect, 
three  items  should  be  included  in  the  analysis.  First,  the  elements  must  be 
curved  elements  for  curved  surfaces,  therefore,  the  displacements  should  be 
functions  of  the  curvilinear  coordinates.  Second,  the  displacements  must  be 
continuous  at  the  edge  of  the  elements.  Third,  the  displacements  and  their 
first  derivatives  should  be  continuous  within  the  domain  of  the  element. 

These  measures  can  be  met  if  the  displacements  are  divided  into  two  por¬ 
tions,  the  displacement  of  the  reference  surface  and  the  displacement  of  a 
director  which  is  initially  normal  to  the  undeformed  reference  surface.  By 
using  a  linear  shape  function  for  the  reference  surface  displacements  and 
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letting  normals  deform,  the  displacements  at  the  edge  of  the  elements  will 
be  continuous. 


To  handle  large  displacements  either  a  Lagrangian  or  an  Hulerian  formula¬ 
tion  must  be  used.  Since  composite  material  properties  are  anisotropic,  the 
calculation  of  material  properties  for  an  Eulerian  formulation  would  be 
complex.  The  properties  of  the  material  are  known  and  established  in  the 
undeformed  configuration.  Therefore,  a  Lagrangian  formulation  will  be  carried 


To  develop  the  strains  and  displacements,  one  must  establish  the  metrics 
in  the  undeformed  and  deformed  systems.  Let  the  position  of  a  point  in  the 

A 

undeformed  shell  be  described  by,  R  (fig  3.1),  where, 

*  *  r  (t)  +  £*  J  (f*)  (3-1) 


(3.1) 


The  position  vector  of  a  point  on  the  reference  surface  is  r  and  *t  is  the 
director  along  the  normal  to  thereference  surface  at  position  r  .  The  basis 


vectors  are  then  given  by. 


ft*  5 
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(3.2) 


where  the  comma  stands  for  covariant  differentiation.  The  metric  of  the 


undeformed  surface,  g^,  is, 


9y  ^  /•  * 

The  position  of  the  same  point  in  the  deformed  shell  is  R  *,  where 

it  — 

n  -  R  +  U 


(3.3) 


(3.4) 


where  U  describes  the  displacement  field.  The  metric  of  the  deformed  body, 
g^j*,  is  then 

q*  =  ft*. 


t 


FIGURE  3.1:  Deformation  Terms  and  Geometry 


(3.5) 
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The  Lagrangian  strain  is  then  defined  as  the  difference  between  the 


metrics. 


(rf- -fr) 


(3.6) 


The  displacement  field  in  its  two  parts  can  be  represented  by, 
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where  U  is  associated  with  the  deformation  of  the  reference  surface  and 
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U  is  associated  with  the  deformation  of  the  director.  Using  equations 
(3.1)  and  (3.7)  and  letting  C«s  <SW  >  where  the  are  the  Lame’  constants 
of  the  surface,  the  strain-displacement  relations  become, 


+J-(UsQJ*(QjsO‘  U,J  (3'8) 
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Now  assuming  the  curvilinear  coordinates  are  the  principal  directions  and  that 
the  transverse  displacement  is  independent  of*  ,  the  strain-displacement  rela¬ 
tions  become,  (see  Appendix  I) 
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(3.11) 
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where 

/ 

R-|=  principal  radius  of  curvature  for  if  direction, 

P  *** 

R2=  principal  radius  of  curvature  for  £  direction, 
a^=  Ae-j  ,  where  e-j  is  unit  vector  in  Jf  direction, 
a*2=  Be2,  where  e2  is  unit  vector  in  jf*  direction, 

3(f)' 0  (?)  e,*t(?)e,*w  (?) 3 

Ci(?)-ti(  ?,?)Z +  *(?,?)£ 

f- 1 

These  Lagrangian  strain  components  then  are  equivalent  to  the  Green 
strain  tensor.  The  stress-strain  law  must  relate  a  Lagrangian  stress  to 
these  Lagrangian  strains.  Since  it  is  Lagrangian  and  symmetric,  the  second 
Piola-Kirchoff  stress  is  used  in  this  work.  This  symmetry  aids  in  estab¬ 
lishing  the  constitutive  law.  The  disadvantage  in  choosing  this  stress  is 
that  it  does  not  have  a  physical  meaning  for  large  displacements.  The  con¬ 


stitutive  law  is, 

S=c  r 


(3.15) 
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represents 


where  S^-  represents  the  second  Piola-Kirchoff  stress  and 
a  fourth  order  constitutive  tensor. 

Since  the  stress  tensor  is  symmetric,  its  order  can  be  reduced.  This 
reduction  is  represented  in  the  conventional  notation  by 

Sll=  S1 

S22  =  s2 

S33  =  s3  (3.16) 

S23  =  S32  =  s4 
S13  =  S31  =  s5 
S12  =  S21  =  s6 


The-order  of  the  constitutive  tensor  and  the  strain  tensor  can  be  similarly 
reduced.  In  addition,  the  components  of  the  strain  tensor  need  to  be  changed 
from  tensorial  to  physical  components,  so  that  the  experimentally  determined 
material  constants  can  be  used.  This  change  is  represented  by 


K 

ixk — 

K  8( /+*//?.) 
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0  _4a _ - 


(3.17) 


The  new  constitutive  law  then  becomes 


The  constitutive  law  for  a  composite  material  takes  on  a  specific  form. 


For  a  material  in  which  the  principal 
ric  axes  coincide,  this  law  in  matrix 


'o„ 

0,i 

0.1 

0 

o 

o 

ott 

O 

0 

o 

D., 

Ot3 

o33 

0 

o 

d 

O 

o 

o 

o 

0 

O 

o 

O 

0 

Orr 

o 

0 

O 

o 

0 

o 

i  case  where 

physical  and 

~o„ 

Qj 

0 

0 

0,o  " 

Oix 

Ox  3 

O 

0 

Dxo 

D,, 

0x3 

0„ 

o 

o 

Ox 

o 

c 

o 

ovv 

Co 

o 

o 

o 

Oys  Off 

O 

► 

0 

O 

Dec. 

material  axes  and  the  principal  geomet- 
form  is 

(3.19) 

geometric  axes  do  not  coincide  becomes, 


(3.20) 


The  physical  components  of  strain  can  also  be  represented  in  matrix  form, 


(3.21) 


The  L  -  operator  matrix  can  be  broken  down  into  1  inear, Lq , and  nonlinear, 
L-j  ,  parts  where 

[LHt.Mr.WM  (3.22 
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The  displacements  (J,  V,  W,  U,  and  V  do  not  explicitly  appear  in  the  LQ 

A- 

matrix.  The  displacements  U,  V,  and  W  appear  to  the  first  power  in  the 

^  A  A  A 

L-j  matrix,  which  is  independent  of  U  and  V.  The  displacements  U  and  V 
appear  to  the  first  power  in  the  L-j  matrix  which  is  independent  of  U,  V 
and  W.  Equation  (3.20)  then  becomes 

S,htDUL.*r.+LliU}  (3.26) 

FINITE  ELEMENT  FORMULATION 

The  theory  and  philosophy  behind  finite  element  analysis  are  well  estab¬ 
lished  and  can  be  found  in  several  sources  and  will  not  be  reproduced 
here  (33  -  35).  The  finite  element  will  be  formulated  using  a  variational  method 
and  a  displacement  model.  First  the  nodal  degrees  of  freedom  and  the  type 


of  interpolation  function  for  these  degrees  of  freedom  must  be  selected. 

The  interpolation  function  must  meet  two  criteria.  It  can  be  one  degree 
less  than  the  order  of  the  governing  differential  equation,  and  it  must  be 
smooth  within  the  domain  of  the  element.  From  Equation  (3.21),  (3.23), 

(3.24),  and  (3.25),  the  order  of  the  differential  equations  is  two;  there¬ 
fore  a  first  order  interpolation  is  sufficient.  For  the  function  to  be  suffi¬ 
ciently  smooth  within  the  domain  of  the  element,  at  least  all  first  deriva¬ 
tives  must  be  continuous.  An  additional  consideration  in  the  case  of  laminated 
composites  is  that  the  displacement  function  be  able  to  change  from  layer  to 
layer.  This  means  that  in  order  for  the  first  derivatives  to  be  continuous 
through  the  thickness,  they  must  be  at  least  linear  within  each  layer. 

A*  A* 

To  meet  these  criteria,  U,  V,  and  W  must  be  at  least  linear  with  respect 

-*  A  A 

to  the  *.  and  X  coordinates,  and  U  and  V  must  be  at  least  linear  with  respect 
to  i  ,X  ,  and  Z. 

24 


Since  the  thrust  of  this  work  is  the  nonlinear  analysis  of  plates  and 
shells  under  varying  loading  and  boundary  conditions,  the  most  convenient 
element  to  use  is  the  triangular  element.  The  use  of  triangular  elements 
allows  better  matching  of  irregular  geometric  boundaries.  The  node  number 
ing  and  geometric  arrangements  are  as  shown  in  Figure  3.2. 

Then  the  interpolation  function  for  the  displacements  associated  with 


the  deformation  of  the  reference  surface  is, 

u-  N,  u,  +  Nx  ut  +N3  u3  ) 
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(3.27) 


(3.28) 


(3.29) 
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(3.31) 


(3.33) 


£-*',)?  (3-32) 
t/y  ax-  r;  O  2  <3-33) 

where  is  the  coordinate  in  the  i  direction  at  node  i. 

To  meet  the  smoothness  and  continuity  conditions  at  the  lamina  inter¬ 
faces,  the  displacements  are  not  only  forced  to  be  equal,  but  also  the  Z- 
derivatives  of  the  displacements  are  forced  to  be  equal.  Defining  the  Z- 
coordinate  of  the  bottom  of  the  ith  lamina  as  Z ■  -j  (Fioure  3.2)  and  the 
Z-coordinate  of  the  top  as  Z-,  the  derivatives  at  the  interfaces  can  be 


defined  as. 


.  d 

^  \  ch 
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(3.34) 

(3.35) 
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Using  equations  (3.34)  through  (3.37),  and  using  a  linear  interpolation 
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(3.38) 

A-  -  A.-I  . 

(3.39) 

To  obtain  the  interpolation  function  for  the  displacements,  a  Taylor  series 
is  expanded  about  the  bottom  of  the  layer  and  truncated  after  three  terms. 

\  '  ...  y  , 
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equations  3.38  and  3.39,  equations  3.40  and  3.41  be 
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where  the  C's  are  constants  which  allow  matching  of  displacements  at  layer 

4  A 

boundaries.  The  constants  can  be  evaluated  by  forcing  U  and  V  to  be  zero  at 

A  4 

the  reference  surface  and  U  and  V  to  be  equal  at  interfaces.  For  this  analy 
sis,  the  reference  surface  will  always  be  the  interior  surface  which  is 
determined  by  the  direction  of  the  surface  normal.  Using  these  criteria, 
equations  (3.42)  and  (3.43)  become, 


jjjL*'  ,d  T*  +  <t)T* 
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where, 
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(3.44) 

(3.45) 
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(3.46) 


where  tj  is  the  thickness  of  the  layer. 


This  manipulation  accounts  for  the  Z-dependence  of  U  and  V.  The  t,  and 
^dependence  is  carried  in  ^  and  9  ,  that  is 

x  <t>-  ^(i*-) 
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To  meet  the  smoothness  and  continuity  requirements,  linear  interpolation 
functions  will  be  used  for  $  and  9.  The  functions  are 

$>-~k (Mi.  <*,  '  NzJ>z  +  Max.  $3  )  (3-48) 
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(3.49) 


The  nodal  degrees  of  freedom  are  the  three  displacements  U,  V,  and  W 


at  each  node  and  the  rotations,  .  $  and  i  Q  ,  at  the  lamina  interfaces  on  the 
director  emanating  from  each  node.  The  element  degree  of  freedom  vector  is 


wktiL  KJjM-ixr  t>f  UujtrS 


(3.50) 


Using  equations  (3.28)  through  (3.30),  (3.44),  (3.45),  (3.48),  and  (3.49), 
the  displacements  in  terms  of  the  nodal  degrees  of  freedom  can  be  represented 
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Equation  (3.21)  then  becomes, 

^e5  =  [L0ft,+L,][w+Nj  fcJ, 

Equation  (3.26)  then  becomes, 


(3.54) 


c^D-'v^jvr0  •  uu  ■  jc  ^  (3.55) 

The  variational  principle  that  will  be  used  to  derive  the  element  stiff¬ 


ness  matrices  is  Hamilton's  principle. 
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^  -C  (^“  dd-  +  ($^££‘0  (3.56) 

where  ^*is  the  kinetic  energy,  Vis  the  potential  energy,  and^K  is  the 
work  done  by  external  forces.  This  analysis  will  be  a  static  analysis; 


therefore,  the  kinetic  energy  is  zero  and  there  is  no  time  dependence.  Equation 
(3.56)  becomes,  -  ~  0  Of 

‘S'W*  (3.57 ) 

Assuming  only  conservative  forces  are  applied  and  letting  "b  or  b1  be  the 

— i  i 

body  forces  and  f  or  f  be  the  surface  forces, 

K.*  £  $}  fal7  ffl d5  0.58) 


(3.58) 


then  for  an  element 


’Md.V+  (3.59) 

From  equation  (3.51), 

$  5 i^rL^  (3. 60 


A  1  T 


(3.6C) 


Equation  (3.59)  can  then  be  written  as, 

rne  A)  r&  rtf&sh  .61 ) 


Now  let, 


(3.62) 


Combining  equations  (3.61)  and  (3.62)  the  variation  of  the  external  work  on 


one  element  is. 


The  potential  energy  for  an  element  can  be  represented  by 


v-iiuruuo 


The  variation  of  the  potential  energy  is 


IcTH^W 


(3.63) 


(3.64) 


(3.65) 
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Using  equations  (3.54)  and  (3.55),  equations  (3.65)  becomes 

+  (3.66) 

The  terms  within  the  integration  sign  then  become  the  element  stiffness 
matrix.  The  stiffness  matrix  can  be  divided  into  linear  and  nonlinear  parts. 


such  that 


where 
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(3.67) 


(3.68) 


(3.69) 


Using  equations  (3.57),  (3.63),  and  (3.67),  the  basic  variational  equation  is 

f]  .  5  ^riX  -  Kwf  J  1^3  (3.70) 

Since  the  variations  are  arbitrary  and  independent,  equation  (3.70)  becomes 

£af$  (3.71) 

Since  the  functions  chosen  to  represent  the  degrees  of  freedom  are  continuous 
at  element  boundary  interfaces  ,  the  equations  for  each  element  can  be  added 


together  to  obtain  a  global  set  of  equations, 

fol‘  (V  Ul 


(3.72) 


SOLUTION  PROCEDURE 

The  next  step  in  the  analysis  is  to  solve  the  set  of  nonlinear 
equations  presented  by  equation  (3.72).  The  Newton-Raphson  method  will 
be  used.  Although  derivations  of  this  method  are  available,  it  is  presented 
here  for  clarity. 
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To  solve  equation  (3.72)  an  initial  estimate  of  thevector  U l  is 


required.  Switching  to  tensor  notation,  let  this  initial  estimate  be  a^. 

J 


Further  let 


Kl  -  K 


and 


=  r... 


(3.73) 


NL 


TJ 


.  o 

then  the  residual  force  t  i,  representing  the  degree  to  which  body  force 
equilibrium  is  not  satisfied,  is  given  by: 


(3.74) 


To  obtain  a  better  estimate  of  a-,  say  a j ,  the  residual  force  vector 


is  expressed  in  a  truncated  Taylor  series  expanded  about  the  estimate  a°. 

J 


The  Taylor  series  is  truncated  after  the  linear  term.  The  expression  is: 

(3.75) 

.  -*  IA-  v.  /  ' f  . 

i.  ? 

Now  it  is  assumed  that  the  new  residual  force  vector  will  be  the  null  vector 
Then  equation  (3.75)  becomes 


-  c  r 


where  (  f^Je  '  f  ^ 


(3.76) 


Using  equation  (3.74),  an  expression  for  is  obtained 

^ f  ^  *-V  ^ 


Or  gathering  terms 


(3.77) 

(3.78) 


Since  K..  is  independent  of  the  a-'s,K..  .  =  0.  Equation  (3.78)  becomes 

•  In  J  *  )  J 

(3.79) 

All  terms  on  the  right  hand  side  of  this  equation  are  known  except  for  the 
third  order  tensor  K-.  . . 

1  K  >  J 

This  tensor  is  evaluated  in  one  of  two  basic  ways.  First,  the  third 
order  tensor  can  be  evaluated  and  the  indicated  tensor  multiplication  carried 
out.  A  second  alternative  is  explained  by  Zienkiewicz  (33).  He  condenses 
Kik  ^a^  into  a  second  order  tensor  called  the  initial  stress  or  geometric 
matrix.  Both  of  these  methods  require  additional,  extensive  calculations  on 
the  elemental  level  for  a  classical  three  dimensional  analysis.  For  an  element 
such  as  described  here,  the  calculations  would  be  even  more  extensive  due  to 
the  coupling  between  the  translational  degrees  of  freedom  and  the  larger  number 
of  rotational  degrees  of  freedom. 

To  avoid  these  additional  calculations,  equation  (3.79)  is  replaced  by 


the  approximation 


(3.80) 


The  inversion  process  is  then  carried  out  as  shown  in  equation  (3.76).  The 
approximation  given  by  equation  (3.80)  will  reduce  the  rate  of  convergence 


(63),  but  this  should  be  offset  by  fewer  calculations  required  at  each  iter¬ 
ation. 

This  process  is  continued  until  either  the  error  as  indicated  by  the 
size  of  the  residual  force  vector  or  the  difference  between  successive  approx¬ 
imations  to  the  degree  of  freedom  vector  is  small.  There  are  various  measures 
that  can  be  used  to  define  the  size  of  these  vectors  and  hence  convergence  (83). 
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The  method  which  is  used  in  this  study  varies  per  degree  of  freedom.  If  the 
degree  of  freedom  is  one  at  which  an  external  force  is  applied,  convergence 
occurs  if  the  size  of  the  residual  force  is  some  set  proportion  less  than  the 
applied  external  force.  If  the  degree  of  freedom  does  not  have  an  externally 
applied  force,  the  size  of  the  residual  force  is  compared  to  an  absolute  value. 
Convergence  occurs  if  the  absolute  value  of  the  residual  force  is  less  than 
the  set  value.  Global  convergence  occurs  if  the  residual  forces  for  all  degrees 
of  freedom  are  within  their  respective  convergence  criteria. 

After  convergence  occurs  at  one  load  level,  the  load  is  proportionally 
increased,  and  the  basic  solution  process  is  repeated  until  a  convergent 
solution  is  reached.  In  previous  nonlinear  analyses,  load  incrementation  is 
continued  until  the  collapse  load  is  reached,  that  is  until  the  stiffness  matrix 
becomes  singular.  The  methods  which  analyze  behavior  in  the  post  collapse  phase 
normally  switch  to  displacement  incrementation. 

The  current  stiffness  parameter  has  been  used  previously  but  a  detailed 
derivation  is  not  available  in  the  literature.  Such  a  derivation  is  presented 
in  Appendix  3  . 

The  basic  formulas  for  the  current  stiffness  parameter  are  as  follows. 

The  stiffness  parameter,  S  ,  is  used  when  the  loading  is  proportional ,  that 
is  when  some  load  fq£  is  defined  by: 

(3.8D 

where  p  is  the  proportional i ty  factor  and  ^ is  a  reference  load,  usually 
the  initial  load.  The  current  stiffness  parameter  is  defined  by. 
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where  K\  is  the  degree  of  freedom  vector  at  the  reference  load  and  i 
is  the  degree  of  freedom  vector  at  the  corresponding  load  fq  |  .  The  partial 
derivatives  are  then  approximated  using  finite  differences.  Equation  (3.82) 
becomes , 


(3.83) 


If  the  load  steps,  p,  are  all  equal  then 

5p  Z  A 

f  A 


(3.84) 


The  sign  of  Sp  indicates  the  stability  condition  of  the  structure 


according  to  the  following  relations. 


Sp>  0  stable  equilibrium 


Sp  =  o  neutral  equilibrium 


(3.85) 


Sp<  0  unstable  equilibrium 


To  apply  equations  (3.85),  it  is  necessary  to  project  the  value  of  Sp.  This 


is  done  with  a  two  term  Taylor  Series  in  which  the  first  derivative  is 
approximated  by  a  first  order  finite  difference  formula 


oy  a  first  order  finite  difference  formula 


(3.86) 


If  the  projected  value  of  Sp  is  negative  and  the  solution  does  not  converge, 


the  point  of  collapse  has  been  reached. 

The  overall  solution  procedure  for  a  nonlinear  problem  is  then: 

1.  Solve  the  linear  problem  and  use  the  linear  solution  as  a  first 
approximation  to  the  solution  of  the  nonlinear  problem. 

2.  Solve  the  nonlinear  problem. 

3.  If  convergence  occurs,  continue;  if  non-convergence  occurs,  stop 
the  solution. 

4.  Calculate  the  current  stiffness  parameter  and  use  eqn.  (3.85)  to 
determine  the  current  state.  Stop  at  neutral  equilibrium. 
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5.  Increment  the  load  by  the  set  proportion  and  recalculate  the  load  vector. 

6.  Return  to  step  2  and  repeat. 
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4.  RESULTS 


In  this  chapter,  three  areas  will  be  covered  which  demonstrate  the  con¬ 
vergence  characteristics  of  the  element,  the  accuracy  of  the  element,  and  the 
applicability  of  the  element  to  nonlinear  laminated  plate  problems. 

The  effects  on  convergence  of  different  conditions  will  be  investigated. 

The  number  of  elements  in  the  reference  surface  wiTT  be  varied  to  examine  the  influ¬ 
ence  of  element  size.  Since  this  element  includes  transverse  effects ,  .the 
number  of  layers  in  the  thickness  direction  of  the  plate  will  also  be  varied. 

Last,  the  effect  of  different  materials  and  plate  geometries  will  be  inves¬ 
tigated.  Five  specific  linear  cases  are  inspected  to  cover  a  range  of  dif¬ 
ferent  materials  and  geometries. 

The  accuracy  of  the  element  is  determined  by  comparing  computer  results 
against  other  solution  results.  The  displacement  accuracy  of  the  element  for 
linear  problems  can  be  inferred  from  the  convergence  study.  Stress  accuracy 
for  linear  problems  is  examined  by  comparing  with  results  obtained  for  three 
semi-infinite,  laminated,  composite  plates.  The  accuracy  of  the  nonlinear 
solution  is  studied  by  comparing  with  the  previously  calculated  results  for 
a  square  orthotropic  plate. 

To  demonstrate  the  range  and  applicability  of  the  element,  two  additional 
nonlinear  problems  will  be  studied.  First  the  change  in  stresses  for  an  ortho¬ 
tropic  semi-infinite  plate  undergoing  large  displacements  will  be  studied. 

Second,  the  method  will  be  used  to  determine  the  collapse  load  for  a  square 
orthotropic  plate  with  a  centrally  located  hole. 


CONVERGENCE 


The  first  case  that  will  be  examined  for  convergence  is  a  semi - i nf i ni te 
isotropic  plate  under  uniform  loading.  The  ratio  of  the  thickness  Ij  the 
width  for  this  plate  is  .02;  therefore  it  can  be  classified  as  a  thin  plate. 

An  illustration  of  the  structure  and  a  typical  finite-element  mesh  are  in 
figure  4.1.  The  graph  showina  the  effects  of  increasing  the  number  of 
elements  or  the  numberof  layers  is  figure  4.2.  Represented  on  the  vertical 
axis  is  the  ratio  of  the  displacement  in  the  z  d i rect ion,w(at  the  midpoint), 
found  using  this  element  and  the  classic  value  of  w  (84).  Plotted  along 
the  horizontal  axis  an  the  divisions  along  the  side  of  the  plate. 

Increasing  the  number  of  divisions  in  a  consistent  manner  causes  the 
results  to  approach  the  classic  result.  Also  increasing  the  number  of  layers 
improves  the  accuracy  of  the  solution.  For  the  plate  with  50  side  divisions 
and  100  elements,  increasing  the  layers  from  one  to  two  improves  the  accuracy 
by  122;  increasing  from  two  to  four  layers  improves  accuracy  by  2%,  and  in¬ 
creasing  from  four  to  six  layers  improves  accuracy  by  .2%.  Any  larger 
number  of  divisions  does  not  noticeably  improve  accuracy.  The  highest  accur¬ 
acy  achieved  was  862  of  the  classical  result  for  a  plate  with  50  side  divisions 
and  8  layers.  If  the  finite  element  solution  continued  to  approach  the  classi¬ 
cal  solution  at  the  same  rate,  220  elements  would  be  required  to  achieve  1002 
accuracy.  The  loading  for  all  cases  is  found  using  eqn.  (3.62) 

The  second  case  to  be  investigated  is  a  square,  clamped,  isotropic  plate 
under  uniform  loading.  A  drawing  of  the  plate  and  a  typical  element  mesh  are 
shown  in  figure  4.3.  The  graph  of  the  results  is  figure  4.4.  The  axes  on  the 
graph  are  the  same  as  those  in  figure  4.2.  The  convergence  rate  relative  to 
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FIGURE  A. 3:  Square  Plate  Geometry  and 
Finite  Element  Model 


FIGURE  4.4:  Convergence  Graph  for  Square  Isotropic  Plate. 
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the  number  of  side  divisions  is  the  same  for  the  semi-infinite  and  the  square 
plate  as  the  number  of  divisions  is  increased.  However  the  two  dimensional 
nature  of  the  square  plate  may  lead  to  a  requirement  for  more  elements.  In 
order  to  obtain  the  classical  results{85),  110  slide  divisions  dnd  approximately 
24,000  elements  v/oul  d  be  required.  Increasing  the  number  of  layers  for  the  square 
plate  did  not  have  a  large  effect  because  the  mesh  was  large. 

The  third  case  is  a  semi-infinite  orthotropic  plate  under  transverse  sinu¬ 
soidal  loading,  q,  where 

q  =  q  sin  jL>L  (4.1) 

u  cu 

The  plate  set-up  and  the  mesh  arrangement  are  the  same  as  shown  in  Figure  4.1. 

The  fibers  run  in  the  x-direction.  The  layers  are  all  of  the  same  thickness 
and  the  ratio  of  plate  thickness  to  width  is  .02.  The  material  properties  are 
given  in  Table  4.1. 


TABLE  4.1  COMPOSITE  MATERIAL  PROPERTIES 


El/  Et  =  25  Glt/  Et  -  .5 


GTT/  ET  "  -2  VLT  VTT  ~  -25 

L  signifies  the  direction  parallel  to 
the  fibers. 


T  signifies  the  direction  perpendicular  to 
the  fibers. 


The  axes  for  the  results  graph.  Figure  4.5,  are  labeled  the  same  as  they  were 
in  the  previous  two  examples.  The  results  are  compared  with  the  analytical 
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results  obtained  by  Pagano  (3). 

The  finite  element  solution  for  the  orthotropic  plate  converges  twice  as 
fast  as  the  solution  for  the  isotropic  plate.  For  the  plate  with  20  side 
divisions,  going  from  one  layer  to  two  layers,  increases  accuracy  by  14.1%, 
going  from  two  layers  to  four  layers  increases  accuracy  by  1.4%,  going  from 
four  to  six  layers  increases  accuracy  by  .1%,  and  going  from  six  to  eight 
layers  increases  accuracy  .01%.  The  most  accurate  solution  is  obtained  with 
20  side  divisions  and  eight  layers  when  the  finite  element  result  is  99%  of 
the  analytical  result. 

The  fourth  problem  that  is  examined  for  convergence  is  a  thick  semi¬ 
infinite  orthotropic  plate.  This  problem  is  similar  to  the  previous  one  except 
that  the  ratio  of  thickness  to  width  is  .25.  As  can  be  seen  from  Figure  4.6, 
the  finite  element  results  converge  quickly  to  the  classical  results  (3). 

With  four  layers  and  four  side  divisions,  the  finite  element  solution  is  equal 
to  the  classical  result  for  the  w  displacement.  Increasing  the  number  of  side 
divisions  past  four  does  not  affect  the  accuracy.  With  four  side  divisions, 
going  from  one  layer  to  two  layers  improves  accuracy  by  13.9%  and  going  from 
two  layers  to  four  layers  improves  accuracy  by  .5%. 

The  last  convergence  example  is  for  a  square,  laminated  composite  plate 
under  transverse  sinusoidal  loading,  q, where 

q  =  q  sin  ^5—  sin  HX  (4.2) 

M  Mo  a.  cu 

The  basic  set-up  of  the  plate  and  the  finite  element  mesh  are  similar  to  that 
shown  in  Figure  4.3.  The  boundary  conditions  for  this  plate  are  simple  supports 
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Orthotropic  Plate 


46 


reproduced  on  usafa  copier 


Along  the  edges  the  transverse  displacement,  s,  is  fixed  and  the  normal  to  the 
reference  surface  at  the  edges  is  not  allowed  to  rotate  perpendicular  to  the 
edge.  The  plate  is  a  symmetric  (0,90)  laminate  of  graphite  epoxy;  the  material 
properties  are  shown  in  Table  4.2.  The  thickness  ratios  vary  from  .02  to  .2. 


As  can  be  seen  in  Figure  4.7,  the  finite  element  solution  for  the  thick 
plates,  where  the  ratio  of  thickness  to  width  is  .2  and  .15,  approaches  the 
shear  deformation  theory  (SDT)  solution  (23)  quickly.  With  five  side  divi¬ 
sions  the  finite  element  solution  is  99%  of  the  SDT  solution.  As  the  plate 
becomes  thinner,  the  rate  of  convergence  decreases.  For  a  plate  with  a  thick 
ness  ratio  of  .02,  a  convergent  solution  is  rot  obtained.  Although  the  rate 
of  convergence  does  appear  to  be  steady. 

Through  these  five  examples,  convergence  has  been  studied.  The  examples 
have  indicated  the  effect  of  the  number  of  elements  and  the  number  of  layers 
on  convergence.  Plates  with  simple  and  clamped  supports  have  been  studied. 
The  materials  have  ranged  from  isotropic  to  laminated  orthotropic.  Now  that 
convergence  has  been  studied,  the  accuracy  of  the  element  for  linear  and  non¬ 
linear  problems  must  be  examined,  -although  the  validity  of  the  method  can  no 
longer  be  questioned.  With  this  method  it  is  possible  to  include  transverse 
shear  effects  and  obtain  usable  results. 
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IDE  DIVISIONS 


ACCURACY 


Although  the  transverse  displacement  accuracy  for  some  linear  cases 
was  examined  in  the  convergence  section,  the  in-plane  displacement  and  stress 
accuracy  for  linear  problems  was  not  investigated.  Three  further  linear  cases 
will  be  studied.  All  three  cases  are  semi-infinite,  anisotropic  plates  under 
sinusoidal  transverse  loading,  equation  4.1.  The  geometry  of  the  problems  are 
shown  in  Figure  4.1.  All  cases  have  the  edges  with  simple  supports.  The 
material  properties  are  given  in  Table  4.1.  Case  1  is  an  orthotropic  plate 
with  the  fibers  aligned  in  the  x  direction.  Case  2  is  a  two  layer,  unsymmetric 
composite  plate  with  the  fibers  in  the  bottom  layer  parallel  to  the  x  axis  and 
the  fibers  in  the  top  layer  parallel  to  the  y  axis.  Case  3  is  a  three  layer 
symmetrically  laminated  composite  plate.  In  the  upper  and  lower  layers,  the 
fibers  are  aligned  in  the  x  direction.  In  the  middle  layer,  the  fibers  are 
aligned  with  the  y  axis.  The  results  of  this  analysis  are  compared  against 
those  of  Pagano  obtained  using  elasticity  (3).  For  all  of  these  cases,  the 
following  normalized  quantities  will  be  used 


-  ..  % *(o,i) 
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For  case  1,  the  accuracy  of  the  w  prediction  over  a  range  of  thick¬ 
ness  ratios  is  examined.  As  can  be  seen  in  Figure  4.8,  the  predictions 
are  accurate  over  a  wide  range  of  S  values.  The  largest  error  is  3%  for 
S=50.  As  can  be  seen  in  Figure  4.9,  the  stresses  are  also  predicted 
accurately  for  a  thick  plate  with  S=4.  The  largest  errors  occur  at  the 
top  and  bottom  of  the  plate.  In  Figure  4.10,  the  results  for  a  slightly 
thinner  plate  are  shown.  The  finite  element  results  match  the  elasti¬ 
city  results  extremely  well.  The  last  graph  for  case  1  is  shown  in 
Figure  4.11.  In  this  graph  the  values  of*E  for  the  elasticity  solu¬ 
tion  and  the  finite  element  solution  are  shown.  The  results  using  this 
finite  element  almost  exactly  match  the  elasticity  results.: 

Since  case  2  represents  an  unsy;,. metrical  laminate,  it  represents  a 
critical  test  for  a  proposed  finite  element.  Comparing  values  for  w  in 
Figure  4.12,  the  finite  element  and  elasticity  results  are  the  same  over 
a  range  of  thickness  ratios.  In  Figure  4.13,  the  values  foro^  are 
compared.  Again,  the  two  results  are  essentially  the  same.  The  next 
area  examined  is  the  transverse  shear  stress,  T  .  Looking  at  Figure 
4.14,  again  the  two  methods  give  the  same  results.  The  largest  differ¬ 
ence  occurs  at  the  bottom  of  the  plate.  The  last  parameter  that  is  in¬ 
vestigated  is  u  for  S=4.  Once  again  looking  at  Figure  4.15,  finite 
element  and  elasticity  results  are  the  same. 

The  last  linear  case  for  the  semi -infinite  laminated  plate  is  case 
3.  As  can  be  seen  from  Figure  4.16,  the  finite  element  prediction  for 
w  loses  accuracy  as  the  plate  becomes  thinner.  This  lack  of  accuracy  is 
caused  by  the  same  phenomena  that  caused  slow  convergence  for  a  thin 
isotropic  plate.  For  the  thin  isotropic  plate  and  the  thin,  laminated, 
symmetric  plate  the  transverse  shear  strain  should  go  to  zero,  whereas 
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FIGURE  4.10:  Semi -infinite  Plate  Case  1:  <7*  vs.  z,  S  =  10 
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for  thick  plates  and  laminated  unsymmetric  plates,  the  shear  strain  is 
non-zero.  The  lack  of  accuracy  results  from  the  choice  of  the  shape 
functions.  From  equations  (1.57,58)  the  shear  strain  depends  on  the 
planar  derivative  of  w  and  the  z  derivative  of  u  and  v.  Since  w  varies 
linearly  with  respect  to  x  and  y,  equation  (3.29),  the  w  derivative  is 
a  constant  within  the  element.  Since  u  and  v  vary  quadratical ly  with 
z,  equations  (3.42,43),  the  derivative  is  linear.  Therefore,  no  matter 
how  small  the  element  is  made  the  transverse  strains  will  never  go  to 
zero  within  the  element,  and  this  element  as  derived  will  never  be  able 
to  accurately  handle  these  cases.  To  resolve  this  problem,  more  bending 
could  be  introduced  in  the  element;  that  is,  if  the  w  shape  function  is 
made  to  be  quadratic,  the  planar  derivatives  would  be  linear,  and  the 
element  could  resolve  itself  to  zero  shear  strain. 

For  case  3,  the  results  for  <r ^  are  nearly  identical  for  S=4  and 
S=10,  Figures  4.17  and  4.18.  The  results  for  the  shear  stress,  ^xz, 
are  also  nearly  identical.  The  small  variation  is  caused  by  the  differ¬ 
ence  in  the  conditions  which  are  satisfied  at  the  layer  interfaces.  In 
Pagano's  solution ,  ,  V ,  u,  and  w  must  be  equal  at  layer  interfaces. 

In  this  finite  element  solution,  the  compatibility  conditions,  u  and 
u,z>  are  made  equal  at  the  interfaces.  Since  w  is  constant  through  the 
thickness,  the  conditions  which  Pagano  satisfied  are  met  by  the  compati¬ 
bility  conditions,  but  the  compatibility  conditions  are  not  satisfied  by 
making  stresses  equal  at  interfaces. 

To  investigate  the  effect  of  geometry  on  accuracy,  a  square  ortho¬ 
tropic  plate  is  examined.  The  present  element  results  are  compared  with 
an  elasticity  solution  (23)  and  another  finite  element  solution  (45). 

The  loading  is  given  by  equation  4.2;  the  boundary  conditions  are 


FIGURE  4.17:  Semi -infinite  Plate 
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described  in  the  fifth  convergence  example,  and  the  material  properties 
are  given  in  Table  4.2.  As  can  be  seen  in  Figure  4.21,  the  finite 
element  results  with  the  present  element  are  better  than  those  obtained 
by  Pryor  and  Barker.  The  primary  difference  between  the  elements  is  that 
the  normal  remained  straight  in  the  element  used  by  Pryor  and  Barker. 

By  allowing  the  normal  to  both  rotate  and  deform,  the  elasticity  results 
can  be  dupl icated. 

The  last  example  that  will  be  investigated  in  this  section  is  the 
nonlinear  deformation  of  a  square  orthotropic  plate.  Since  few  nonlinear 
analyses  of  orthctropic  plates  have  been  done,  it  is  difficult  to  find 
examples  by  which  accuracy  can  be  judged.  The  finite  element  results 
will  be  compared  against  some  approximate  elasticity  results  obtained  by 
Chia  (77).  The  plate  is  clamped  on  all  edges  and  under  a  uniformly 
distributed  transverse  load.  The  material  properties  are  given  in  Table 
4.2.  The  comparison  is  presented  in  Figure  4.22.  The  finite  element 
results  and  elasticity  results  are  very  close  for  this  example. 

RANGE  AND  APPLICABILITY 

Two  examples  will  be  used  to  illustrate  the  applicability  to  non¬ 
linear  problems.  First,  the  stresses  in  a  semi -infin i te  unsymmet^i c , 
laminated  plate  will  be  examined.  Second  the  collapse  load  for  a  square 
orthotropic  plate  *;ith  a  central  hole  will  be  determined. 

The  semi -i nfi ni te  orthotropic  plate  used  In  this  analysis  is  the 
one  described  in  case  2  in  the  accuracy  porti’oa  of  this  chapter.  The 
loading  and  boundary  conditions  are  also  the  same. 

The  plot  of  normalized  transverse  displacement  versus  the  normal¬ 
ized  load  is  shown  in  Figure  4.23.  As  can  be  seen,  nonlinear  behavior 
is  obtained.  The  load  at  point  one  corresponds  to  a  linear  load.  While 
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FIGURE  4.25:  Finite  Element  <7"  Plot  for  Unsyrmietric , 
Semi-infinite,  Laminated  Plate,  S=40. 
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loads  two  and  three  are  nonlinear  loads. 

Shown  in  Figure  4.24  are  plots  of  the  normalized  transverse  shear 
stress,  <txz>  The  distribution  shape  o^'cxz  does  not  change  as  the  load 
becomes  larger  but  the  magnitudes  do  appear  to  move  in  the  negative 
direction.  The  plots  of  the  normal  stress,?^,  are  shown  in  Figure  4.25. 
The  same  behavior  as  is  observed  for  'C  is  seen  incr. 

The  last  example  ir  this  chapter  is  an  orthotropic  plate  with  a 


central  hole  under  compressive  load.  The  geometry  of  the  plate  and  the 
finite  element  mesh  use  are  shown  in  Figure  4.26. 


The  plot  of  load  versus  displacement  for  a  plate  with  an  S=50  thick¬ 
ness  ratio  is  shown  in  Figure  4.27.  The  load  is  normalized  by  dividing 
it  by  the  buckling  load,  p  ,  for  an  equivalent  orthotropic  plate  with  no 
hole.  The  transverse  displacement  is  normalized  by  dividing  by  the  dis¬ 
placement  that  would  be  experienced  if  the  plate  behavior  were  linear  to 
the  buckling  load.  The  use  of  the  current  stiffness  parameter  allows 
the  program  to  follow  the  behavior  of  the  plate  up  to  the  collapse  load. 
According  to  this  analysis,  the  plate  collapses  when  the  compressive 
load  is  76%  of  the  classical  buckling  load  for  a  perfect  square  plate. 

In  the  Results  chapter,  the  capabilities  of  this  new  element  have 
been  demonstrated.  With  this  element,  it  is  now  possible  to  obtain 


accurate,  convergent  results  for  the  linear  and  non-linear  analysis  of 
plates.  The  boundary  conditions  can  vary  in  the  transverse  direction. 
The  more  specific  accomplishments  of  this  dissertation  are  contained  in 


the  Conclusions  chapter. 
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5.  Conclusions 


The  success  of  this  dissertation  must.be  measured  against  satisfying 
the  original  objective.  The  stated  objective  of  this  work  is  to  formulate 
a  numerical  procedure  for  the  large-displacement  analysis  of  laminated  com¬ 
posite  structures  with  general  loading  and  boundary  conditions.  In  the 
"Introduction"  chapter,  it  was  also  concluded  that  to  accomplish  this  objec¬ 
tive,  nonlinear  geometrical  effects  and  normal  shear  effects  must  be 
accounted  for. 

To  guage  the  success  of  this  research,  the  development  of  the  theory  and 
the  comparisons  with  other  results  are  examined.  Afte^  these  two  areas  are 
examined,  it  is  then  possible  to  make  conclusions  as  to  how  completely  the 
objective  has  been  met.  Last,  in  this  chapter,  I  will  delineate  the  advantage 
of  using  this  method. 

The  "Background"  chapter  and  the  "Theory"  chapter  outline  the  develop¬ 
ment  of  the  method  used  in  this  dissertation.  The  first  choice  that  was  made 
was  between  finite  differences  and  finite  elements  for  the  numerical  procedure. 
The  finite  element  method  was  chosen  because  of  its  greater  flexibility  in 
handling  boundary  conditions.  This  flexibility  was  demonstrated  by  solving 
problems  with  classical  clamped  simple  support,  and  free  edge  boundary  condi¬ 
tions.  In  addition,  a  mixed  type  of  boundary  condition  was  used  where  a  dis¬ 
placement  normal  to  the  edge  but  in  the  plane  of  the  plate  was  free,  but  the 
normal  at  the  edge  did  not  deform.  The  finite  difference  method  would  not 
have  handled  this  variety  of  boundary  conditions  well. 


The  primary  unique  characteristic  of  this  work  was  the  way  in  which 
normal  shear  effects  were  handled.  Previous  methods  have  allowed  the  normal 
to  rotate  but  not  deform  within  each  layer  of  a  laminate.  This  approach 
resulted  in  discontinuities  in  shear  strain  at  lamina  interfaces.  The  finite 
element  presented  here  overcomes  these  discontinuities  by  guaranteeing  not 
only  displacement  but  also  shear  strain  continuity  in  the  thickness  direction. 
The  original  intent  of  this  approach  was  to  better  model  shear  stresses  without 
the  complexity  of  a  hybrid  finite  element.  As  can  be  seen  in  the  results 
section,  shear  stress  calculations  using  this  finite  element  are  almost  the  same 
as  those  obtained  by  classical  elasticity  methods. 

Another  unique  characteristic  of  the  development  presented  in  this  dis¬ 
sertation  is  the  way  in  which  nonlinear  effects  are  included.  The  nonlinear 
effects  are  included  at  the  basic  elastic  equation  level  rather  than  by  extend¬ 
ing  a  generalized  theory.  This  minimizes  the  assumptions  that  are  made 
and  increases  flexibility  in  the  application  of  the  element. 

The  specialization  of  the  element  now  occurs  through  the  geometrical 
definition  of  the  structure.  Using  this  approach,  it  is  easy  to  specialize 
the  element  to  plates  or  shells. 

To  illustrate  how  the  element  is  used  to  analyze  a  specific  structure, 
plate  problems  were  chosen.  Plate  problems  were  chosen  because  there  were  many 
examples  of  laminated  plate  analysis  to  use  in  comparison  and  because  applying 
the  element  to  plates  was  a. logical  first  step  as  opposed  to  jumping  into  the 
analysis  of  shells. 

After  the  element  was  developed  the  next  area  that  was  approached  was 
the  solution  procedure.  A  classical  nonlinear  equation  solution  procedure 
was  used.  There  was  no  analysis  done  on  how  different  solution  procedures 


affect  solution  efficiency;  therefore,  there  is  no  efficiency  analysis 
accomplished . 

The  one  area  af  the  solution  process  that  includes  unique  information 
is  the  use  of  the  current  stiffness  parameter.  This  parameter  is 
used  to  determine  when  the  nonlinear  stiffness  matrix  becomes  singular 
or  otherwise  when  structural  collapse  occurs.  The  parameter  has  been 
presented  before,  but  a  derivation  of  the  parameter  has  not  been  pub¬ 
lished  to  the  best  of  my  knowledge. 

flow  that  the  development  process  has  been  summarized  and  the 
reasoning  behind  the  process  presented,  it  is  necessary  to  examine 
how  solutions  obtained  using  this  element  compare  with  other  numerical 
and  analytical  results.  The  areas  examined  for  comparison  e  con¬ 
vergence,  accuracy,  and  applicability. 

As  is  shown  and  explained  in  the  "Results"  chapter,  the  element 
did  not  converge  well  and  was  not  accurate  for  thin,  isotropic  plates 
and  thin,  laminated,  symmetric  plates.  The  cause  of  this  was  that  the 
shear  strain  would  not  go  to  zero  as  is  physically  expected(42).  Since 
this  element  was  meant  to  handle  cases  where  the  transverse  shear  is 
an  important  factor  this  shortcoming  was  somewhat  expected,  but  was  not 
resolved  in  this  work  other  than  to  suggest  a  possible  solution. 

The  element  converged  well  for’  anisotropic  materials.  The  ele¬ 
ment  converged  well  for  the  range  of  thicknesses  studied,  but  it  con¬ 
verged  faster  for  thicker  plates.  From  examining  the  semi-infinite 
plate  problem,  it  can  be  seen  that  the  element  converged  faster  for 
the  unsymmetric  laminate. 

Increasing  the  number  of  modelling  layers  in  a  lamina  does  in¬ 
crease  the  convergence  rate.  Using  more  than  two  layer-  to  model  a 
lamina  did  not  greatly  .increase  the  convergence  rate.  Kcre  layers 
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stiffness  parameter  facilitates  easier  nonlinear  analysis.  Even  though  these 
observations  are  made  for  plate  analysis,  the  general  character  of  the  develop 
ment  makes  it  possible  to  adapt  this  element  to  other  structural  shapes. 

As  a  final  statement,  I  would  conclude  that  the  objective  of  this  disserta 
tion  has  been  accomplished.  Even  though  the  element  does  not  handle  isotropic 
cases  well,  it  can  be  used  for  the  nonlinear  analysis  of  laminated  composites 
under  a  variety  of  loading  and  boundary  conditions. 
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APPENDIX  I:  Derivation  of  Strain-Displacement  Relations 


The  objective  of  establishing  these  equations  is  to  describe  the  relation¬ 
ship  between  the  displacements  of  a  point  in  a  structure  of  arbitrary, 
curvature  and  thickness  and  the  strain  at  the  same  point.  To  start,  the  location 
of  the  point  must  be  established  relative  to  a  three-dimensional  reference  frame, 
Figure  3.1.  The  location  is  given  by  ,  where 


(K) 


C**l,  2 


(1.1) 


In  this  relation,  rf-  is  a  director  which  is  perpendicular  to  a  reference 
surface  and  passes  through  the  point  located  by  7?  ,  and  r  is  the  location 

of^the  intersection  of  the  director  with  the  reference  surface.  Both  r  and 
d  are  functions  of  only  the  two  surface  coordinates.  The  metric  of  the 
deformed  surface,  ^  ;j  is  then  given  by 


** ' R' 


/,  z,3 


(1.2) 


where  the  comma  indicates  covariant  differentiation. 

— * 

If  the  poi^t  now  undergoes  a  displacement  iU  ,  its  new  location  is 
described  by  R*  ,  where 


4  — ^ 

R  -  R  -h  UL 


The  metric  in  the  deformed  configuration  i: 

54* 

Substituting  1.3  into  1.4  and  differentiating. 


(1.3) 


(1.4) 


***  ■  VV'V"*'  (,-5) 

The  lagrangian  strain  is  then  defined  using  the  metrics  by 
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Using  equations  1.2  and  1.5,  equation  1.6  becomes 


Simp! ifying 


4  *  *  <4/  4 -  K  ■  4  +  ^  f  4  •  4- 


:  *  . 


—i 

Now  let, the  displacement  be  divided  into  two  parts,  Figure  3.1. 
One  part,  U,  ,  which  is  a  function  of  only  the  surface  coordinates  and 
£*  and  is  associated  with  the  deformation  of  the  reference  surface.  The 
other  part,  XL  ,  which  is  a  function  of  all  thrge  coordinates  and  is 
associated  with  the  deformation  of  the  director.  LL  is  then  given  by 


TL'tCt)*-  £c &,£,*) 


4*1,  X 


The  derivative  of  U-  can  then  be  represented  by, 


Uyj 


The  derivatives  of  K  can  then  be  obtained  from  equation  1.1 


v  I  ^  « 

Now  letting  £  cl^  where  CL^  are 

the  Lame 'constants  of  the  reference  surface,  the  strain  equations  1.8 


( 


fast  Hi + §3£,a  )  -0 + v ) * 

*  £„)  + 

4  4  [(i+ftjA )  •  4, ' +  ^ ' ‘<V  ^  } 

_  _i  _>  > 

— *  A  |  A  A 

^33"  ^  "z 

Now  let  (f^  be  unit  vectors  in  the  surface  such  that 


Using  1.17  in  equations  1.14,  1.15,  and  1.16,  the  strains  become 

<^=  z  M  'C&ip'l'  ^,0  )  * 

-cK+^ 


(i 


(i 


GC 


'3 


(i 


Now  it  is  assumed  that  the  reference  surface  coordinates  are  measured 
in  the  principal  directions.  This  assumption  leads  to  two  identities 


:RVihfbPTc 


*  (L-  *  o 


3  v:- 
1  'v 


Further,  the  assumption  is  also  made  that  the  director  only  deforms  and  does 
not  change  length.  Using  this  assumption,  it  is  possible  to  represent  the 
displacements  as  follow: 


aj  /v 


=  tie.  +-  wet 


Cl*  u£,  +  v£z 


Then  the  derivatives  of  these  displacements  seen  in  equation  1 .10 can  be 
wri tten 


r>J  — ^  t*j  . 


“  +  +  V’^  * 


Since  the  basic  vector  and  the  reference  surface  displacements  depend  only 
on  the  surface  coordinate 

\  ^',3  '  C*,3  ^  ^73  "  0  (I 


OL  is  a  unit  vector  which  only  depends  on  surface  coordinates, 


cL  •  cL,  -  0 


Now  the  Lagrangian  strains  are  expressed  as: 

+  £  •  £. ,  r<t,  ts,*  n, ,  K  v*e , -k  1 1 


Mi] 


*e*,i’K(Pi  tfX? cw)  V  5,-  J,(  [i  w(zf+ 

&■=  l°-zC^+K)f  .(d-tOcu, 

*  V  ^  V  it.C^fOCs  <  i *M 


+  <Lj,i-(tXi2.(iL-i- £L)G+ 0)  it, ^cL (CL*  u)u/,  -h e, ^  ^ 

(fl^ute+uO  '  ^  '  * 

*■  s1_-<i2.  Z/v'2+iz/V(3 +  £vgi,1  jjz(V+v)  j 

*  +V' 

+  <V  £z  &  (2-j  *+  3 

z  ?vv  + **&'  ^  V*,  -  ^,Y^i+  o 

*  ( v  ^  *'-*  1  “>X^ 

+  £,-4,z /^+£,/^Xv+W'C'£/kJi+^,z)(^Vw) 

+  ev  ■  £  (Hfu)^^  tlrti,a^e,,r^(iUSX9*9) 

*  &/  et,[H-tcL)0,ii'V,)tiii2.-  d{Hru-) 

->  £/,! ■  ^2, , (CL-m)( v+0)te.,,2--  <t„  ( ^  (% tu)) 
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From  Equation  (1.16) 


Expanding  this  and  using  equation 


fi:28) 


a  i 

+  V 
>b 


(1.16) 


(1.33) 


From  equation  (1.19) 


f3=i  [a., 0,3  +  «i,  sty, »+H  ; 

*  i, .  t,t  a,3  (£+«,)  +e,  •  4,,  S  ^  f">V  ^ 

(jtv)  W,3] 


(1.34) 
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From  Gould  (4)  and  Calcote  (8) 
vectors  can  be  represented  by 


the  derivatives  of  the  basic 


C  -  7- 

V 


e-,s 


*-/'  & 
~7T~  ^ 


p  .  ii£-£ 
ev  ' 


C,-  '  ^ 


cL,t  * 


— S\ 


where  and  R„  are  the  principal  radii  of  curvature  for  the  reference 
surface.  Substituting  equations  (1.36)  through  (1.41)  into  equations 
(1.30)  through  (I.i5)  a  further  expansion  of  the  Lagrangian  strain  is 
obtained,  where  now  q_  -  ft  Cl^  B  ^-2, 

-k  f  [£u,  + 


i  f[  fy(a+i0  +(v  ^lh\y6r^ 

+  [w-  Hv*v)j}  1 


4=  *  K'+r)[<v4)- % WA) 

+/4<"'0'  -^(Ji-tO?iJlo,L+v,$+  ^Lf {1+$!)+  £l  w7 


-g-  ^U,T  l^//  L^v/2.t  v'2/T  (jJL+LL)+  |~  IA/J 

&r  $(X+MJ  fe£ft+0)Jl  "  "-44 


A  \7_  / ~ 


-y-  y^V'3 


0+  <1 


§f  W'3  +  * 

4  •  i  fft(i+t;)  Ur  *&r/z&u)]+£,3l@(,l-*2, )+■$&*$) 


(1.45 
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%  w]+  %C^,-^(u-tH)j*\>0 [(!>,/ u-46 


(1.47 


Equations  (1.42)  through  (1.47)  are  the  tensorial  components  of 
strain.  To  adjust  them  to  the  physical  components  of  strain,  so  that  the 
applicable  physical  constitutive  equations  can  be  used,  they  must  be 
divided  by  the  Lam'e  constants.  In  anticipation  of  using  these  strains 
in  matrix  equations,  I  will  also  rename  and  reorder  the  strains  such  that 


L  _ 

&0+  Vr.) 


z. 


4. 


l^TTTVeJt- 


(1.48 


(1.49 


i 

(1.51) 

(1.52) 

(1.53) 

The  physical  strains  are  then  given  by 
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4'4  '80^2? 
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(1-55) 


.(1.56) 


e,  ^,)  fa,  -  A.  * 

~B  (v  +  C)*  $u/]  *[(%*?  )~  /  (1-54) 

C^oO^J^t, •t  i.) * 4* (p*Ct)  *  fiK  wJ  *  IB'O )  + 

(I‘55) 

Ify/A,  )*+(%/$.> )  *  1  ' ( 1 ' 56 } 

V  k+jfaobtx"**  ( ?*&)] \S0*%) (&* 

4  [( W'-fr  %[k~  i  (*+  ( 1  • 5  7 } 

V  A/ #£%)[%“ 4  (Z*u  )]**(£%)  faj. (6;v/)‘jrfc*0h 

iv)hj(z;ih  (u.til I  (I-58) 

Bd+\ ^j[(AxAJ'  (* * fiy fe ♦  k)~'&~  (S^u)]^ 

/  f  (4/  4,)  *  *-'(r  ♦  /)*  f  i v  Jfc;  «, )  -| u  (?*&)]+ 
&vy*Ww) IUzsvJ*# (**£)+%  w]-* 

[vi/,-  ^  (b  +  u)]  Lwjx-  \  (v*v)l  l 
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APPENDIX  II  -  Derivation  of  Plate  Finite  Element 


From  equation  3.68,  the  linear  part  of  the  element  stiffness  matrix  is 

[kH]s  j l H  +  [  0J[LolLfi+A/]jV  (II. 

■v 

Expanding  this 

[ kcl.  1-  l WTluflDl L u][a/Uv+ 
fjmLTlDnLj[M]jU  + 
l  L^JT[lJ[o] J  v  + 

fyLNlT[LfLo][LJ  [tf  IV  („. 

Using  linear  algebra 

mLj[ol[L»(MWoJUMf  in. 

Then  equation  (II. 2)  becomes 

[Kfb  [tfiiWoULjmjv* 

]mTiaio][umdv  * 

'v 

Uimunoiuji^yfi 

{WTUf  IdIIlJ]  [//]  J v  <"-4 

The  matrices,  LM  and  w  ,  are  given  in  equation  (3.52)  and 


(3.53).  For  a  plate , 


A  =  1 
B  =  1 


2 

i>  *  x 

U  •* 


(II. 5) 


Then 


Looking  first  at  the  integral 


(II. 6) 


+Y  7* 

By  examination  N L-0  will  have  no  dependence  on  X,  Y,  or  Z  with  the 
given  shape  functions.  The  matrix,  0,  has  no  specific  dependence  on  Z, 
but  its  values  can  change  from  layer  to  layer.  Using  these  relations 


lNTLlDleNJ V*  NrL(itDJi JyJ*  <».8) 

The  integral  Lfydyd/  is  the  area  of  the  element  triangle  at 
the  reference  surface  which  will  be  designated  as  4  •  The  elements  of  the 
constitutive  matrix  in  a  layer  k  are  . 

The  integration  over  can  be  accounted  for  by  defining 


(II. 9) 

n  =  number  of  layers 
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The  integral  then  becomes 


(ht  Ll  d 


’)v{Ztl7 JSl„  &)&  3  &  fa 


The  elements  of  A7  with  dimensions  [y+ZC***)  ,  ?  *  l(n*/)2  are 

j  *■  k^)  +  i7y(£,i  *75u^f) 

^4  J-H) 

K]  '  &&  &  '4  &) 

fa;t3J3  *#,&)  *  %1  &) 

^7^„Ty- *#vs&)+ (ffw7f  *Ps?fc) 

^'*hj**  ^•*3/  JVC.  ~  ^  V 

(where  the  t.  are  the  same  as  the  shape  functions  N-j ) 
Looking  now  at  the  integrals 


(II. io: 


are  then 


(II. ID 


(#TiJ  OLt  fhv  +  ($rfiT£o  L,UJb)'  (II12> 

The  sum  of  these  two  integrals  is  a  symmetric  matrix.  The  matrices  A/ 

A. 

andZ.0  do  not  depend  on  Z,  but  the  matrices  O  and  A/  do  have  Z 

dependence.  The  integration  can  be  accounted  for  by  defining  a  new 
consti tui ti ve  matrix  £ 


where 


9r*z  ±fn* 

<c  1  Dy 


ij  =  1,2, 3, 6 

k  =  interface  numbers  and 
layer  number 

n  =  number  of  layers 


Q'-tS.  /ft. 


i.j  =  4,5  (11.13) 


Using  these  definitions 


jfirL7DL0/vjv  =  r* 


(11.14) 


2 

The  elements  of  are 


H  #  *  &  % ) 

X-  fj  £(%  #  *  K$) '  %(K  % '%  %) 


i  +Jf9*3Ktj 


M  ?  £lj.  Jt-r  +  FT  ^  Y‘  ) 

1 '/«  a  (Btr  tv  QrVirJ 


^  ':^3  tn- 


A  5  Q  >/  •  \  * 


where 


/*zj  r  >T  X  ^  ^ 


Using  equation  (II. 3) 


3X7  --  2X?T 


Examining  the  last  part  of  the  linear  stiffness  matrix,  we  now 
integrate  the  expression 


{.LmTiuYcdi  [  l,  ]  [m]Jv 


(11.16) 


The  integration  over  2 
constitutive  matrix  by 


is  accounted  for  by  defining  terms  of  the 


fL"L?  +  7T  ('S£<!'  7?~) 


-f-  £**■ 


/*/),  q+x.*  a*-  9  ^  /jL4t  ^  /  /x ) 


(11.19) 


* 0  1 ... 


si- 


Then  the  total  linear  element  stiffness  matrix  is 

KJ-rrjVi  +  Vi +  Vi  +  4/iJ  (ii-2o) 

From  equation  (3.69)  there  are  three  parts  to  the  nonlinear  stiffness 
matrix 


v/ 
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Looking  first  at  the  portion  composed  of  a  matrix  and  its  transpose 

and  expanding. 

tNTL,TDL,P  ffjrL/DL,/0  i 

NrLj  bL,hl  +  +  Nt LlbL,ti  +/j7ilTt>L,A/?)dV 

+  ^ jNrr DL'N  *-Nrij0 LtHf  NrL,  TbL. NrH^C,  rb Le N 

+  NXrbL.Z+AJTr,rbL./hM,  D 1*$ ^Lf  b L.N^lv 


(11.22) 


Since  the  second  matrix  integral  is  the  transpose  of  the  first,  it 
will  not  be  calculated,  but  the  transpose  of  the  first  integral  will  be 
added  to  the  nonlinear  portion  of  the  element  stiffness  matrix. 


As  stated  in  Chapter  3,  I  have  assumed  that  the  rotation  will  be  small. 


Using  this  assumption,  any  term  containing  ca  can  be  eliminated 
since  this  matrix  contains  the  rotation  terms  and  the  derivatives  of  the 
rotation  terms.  Applying  this  assumption.  Equation  (11.22)  becomes 

ff^LlOLfirNWL,  N  +NTLrbLfi  * 

^  L  ^  A, 


N  L/VL,  N 


(11.23) 


Examining  the  first  part  of  this  i ntegral  ,_/Ajr/.J"DL,  N 5  the  only 


matrix  that  varies  with 


*  is  M  ;  therefore 

fNrL,r  rLj(£ T)<Li ) L, N dycL x 


(11.24) 
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Using  equation  (II. 9) 

/ hi rLj bt  fi Ay-  £$ ftcLytsLx 


where  X  ,  *)£* 


k«,M-  i¥r  &(£*#)**&. 


>y,3 


•v  '-*J 

where  it.-  U. 

n/l 

6U  =  V 

Of  w 

Mow  let  ^  M  -  ^Mr 

The  next  part  of  equation  (II. 23 )  to  be  examined  is 

[N%rbL,hlcLV 

The  integration  over  Z  can  be  accounted  for  by  using  the  elastic 
constants  defined  by  equation  (11.17).  Using  these  definitions,  the 
integral  becomes 

f,N%r[>L,hJd-V,tA  /U 
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i 

! 


Also  let 


i  T 
A) 


i ,j  =  1,2,3 

k,l  =  0,1,  •  •  *  j  n 

(11.30) 


Looking  at  the  next  integral 


(11.31) 


The  Z  integration  is  accounted  for  by  using  the  elastic  constant 
definition,  equation  (11.13).  Then 


{ZrLr,OLlVjV‘ 

'y  o  ' 


/ 


2  A 


(11.32) 
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H„ ,. 4  37^  flf  #  #  **«&'##  <  £  f)  *  # 

T*^  #  *  *ai  $  # '  a:  £  $0 


'  '‘H.  ^+j  *3K  d  y  A  tv  / 


^r  ^  0„'  *  K-  & 


h^^ruij  &(ij£  a  &  £> 

n(i &:  *  is*  fy  ty+Mif  &) 


H„,  ,J6,  „J,  ^7  (t  K  &&+*  K  tf  fy&tf  #) 

&  a  #*if&+i  4;  fy  $*-  -4;  j*  %) 

>W**5  #  *$*«  &K&* 

#  Hr*(v*  %  +Z£r  fy) 


TT^ 


*  ’A""?* '.v  .  *  ,  -  A-  1 


* 


The  last  integral  in  this  series  is 


/>T4.r  Dl,#jV 

Using  (11.13),  accounts  for  the  2  integration. 
Then 

Dl,MAV=  it  A? 

The  terms  of  h  are  given  by 


" )V)  ~  “  <^  u-^ 

j+3/J-i)z 


J US JL\ 

?r*?7 


v  '*• 

where  u,s  c* 
Pt’-  V 
&/-  w 

Further,  let 

,aM=  V)r 


(11.35) 


(11.36) 


(11.37) 


(11.38) 
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The  last  part  of  the  nonlinear  stiffness  matrix  is 


Clfitfifli,  <  Mi,  n,  llvtA/UV 


(II.3C 


Expanding  this  integral  and  eliminating  L,  terms 


fa  amc i  Lv-Z/Oi  j i/~jffirL,  rDL,fi 


H  rOL,i  +  /f'CDl  N  W  TL,TPi:M  J 
M7l,  dl,  d  *(# tC Ol,m) 


(II.4C 


17- 


Si  nee 
(11.40)  becomes 


,  equation 


[UhAfLustfim  UW.  1 LMW* 

j*{  fiTC,  tol,  n  *  'V’T,  tdl,  h*(a vtl,pl,m)t+ 

hit  DL,  v}  JV 


(11.41 


Examining  the  first  term  of  (11.41)  and  using  (II. 9)  to  defin  elastic 
constants 


The  next  term  in  (11.41)  is 


L  fiTLj  DL,  /DjV 


(11.44) 


By  using  equation  (11.13)  and  defining  a  new  matrix  /A  ,  the  result  is 


(11.45) 


where 


K 

it C  #  ft+K 


l 


y 

J.  c^/h  e)Lj_r  f  1 

X 


V*,**1;  v 


To  account  for  the  third  term  in  the  integral 


/  / 

(xx. y/; 


(11.46) 


let 


!S~  /  y  -j- 

/'i*  Vir 


/<Klt 


TheApart  of  (11.41)  to  reduce  is 


ffNTL? o  £,  A IV 


(11-47) 


The  ?  integration  in  this  integral  is  accounted  for  by  using  the 
consti tui tive  matrix  defined  by  (11.17).  The  integral  is  then  represented 
by 

//V  L.OL'AJJl/*  26  >? 

(11.48 


where 


X; 


'?+/  +34  * 3  -l)U  4  7f j  +3 Jr  3(f  ~/)(a  n)***  **  U*^ 


ft 


U  e) 
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AJ  --O  ! ,  “  "  /T 


~  / ,  ~L 
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The  element  nonlinear  stiffness  matrix  is 


[a^7=  it  [rA7  +  (y\  +  nrt  "  *  m  *  *rt  4  ^°m  + 


IT-  -  j3  1 4  if  J(*  1 

y)  -*  *  /?  ■*  yiJ 


(11.49 
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APPENDIX  III:  Derivation  of  Current  Stiffness  Parameter 


The  objective  of  deriving  the  current  stiffness  parameter  is  to  clarify 
the  conditions  which  govern  its  use.  The  current  stiffness  parameter  is  a 
scalar  which  can  be  used  to  predict  collapse  loads  for  continuous  structures. 
By  using  this  scalar,  the  problem  of  examining  the  singularity  of  a  nonlinear 
stiffness  matrix  is  avoided. 


As  a  structure  becomes  unstable,  there  will  start  to  be  an  unlimited  dis¬ 
placement  for  a  small  increase  in  load.  Representing  the  continuous  structure 
as  an  infinite  number  of  degrees  of  freedom,  this  instability  is  represented 
mathematically  as, 


/  =1,2 


( 1 1 1  - 1 ) 


where  a-  is  a  generalized  displacement  and  q.  is  a  generalized  force. 

For  general  loadings,  the  sign  of  cannot  be  determined  prior  to 

collapse.  If  proportional  loading  is  used,  some  conclusions  can  be  drawn. 
Proportional  loading  is  defined  by 


1 

qf  >  P(<V>( 


'=1.2 


(III. 2) 


where  is  the  load  on  the  structure,  p  is  a  proportionality  factor,  and 


Kb 


is  an  initial  load  which  is  less  than  the  collapse  load. 

Even  with  proportional  loading,  the  sign  of  cannot  be  determined, 

but  the  sign  of  will  be  the  same  as  the  sign  of  ~3^l  prior 

to  collapse.  Since  these  terms  are  the  same  sion,  the  ratio  of  these  two  scalars 

<7  cl  l 

will  be  positive^  prior  to  collapse.  Since  ~T$~i  becomes  infinite  as 
collapse  and  ~3(a'f  is  some  finite  number,  collapse  is  defined  by 

cfUJ; 

(III. 3) 


_ =  Q 


The  point  of  instability  is  either  a  local  maximum,  a  local  minimum,  or 
an  inflexion  point  for  the  total  load-di solacement  curve.  If  the  point  of 
instability  is  a  local  maximum  or  minimum,  the  sign  of  the  derivative  ratio 


will  change.  If  the  point  of  instability  is  an  inflexion  point,  the  sign 
of  the  derivative  ratio  will  not  change.  Since  the  derivative  ratio  relates 
the  way  the  structure  is  currently  carrying  a  load  to  the  way  it  originally 
carried  a  load,  it  can  be  called  a  current  stiffness  parameter.  Since  it  is 
good  only  for  proportional  loading,  this  must  be  indicated  in  the  symbol,  S 


df/ 

With  this  scalar,  the  current  stability  status  of  a  structure  can  be 
determined,  that  is 


S  >0  structure  stable  or  past  inflection  point  headed  for 
p  collapse. 

Sp  =  0  collapse  point. 

Sp  0  post  collapse. 

Since  this  structure  is  undergoing  proportional  loading,  S  can  be 
further  simplified.  p 

Using  the  chain  rule 


Then 


V 


d  P 

dp  J(%X 


(II 


e) 

dp  dfy 


Now  multiplying  numerator  and  denominator  by  the  scalar 


#  4 
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t)(°-r)*.  P  c) f>  c) 1  j/c 

sp  =  c> p~~  c>&*  ^a-  ~zy 


where  is  the  Kronectfer  delta. 


From  the  definition  of  proportional  loading. 


Using  this  relation, 


equation  (IV. 7)  becomes 
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'<■*)»  t)p  Jp  J  [fc jf 

f3  +  <)  P  - 

m  7»\^. 

^  ^ tPp  J  (ji')p 

~ZT~  7^  di^r)* 


111 


but  also 


■  &  -  p  % 

c  d&ijZe. 

bD  “  ^  P  * P 


Multiplying  numerator  and  denominator  by  the  scalar  ^  ~Tp 

Jp  dp  J*K  j<lk  Jp  ~Tp 


p 

*.-Cj 

- . - —  r  ■ 

i^A:  Jii  Ji*  > 

A  p  a>jo  J  p 

Jp 

Jp  ^ 

^p'  J 

c 

Je:  r*  r 

Jp  Jp  vn*  ** K 

sp  “ 

since  =  <4cn 

and  j,  k,  m,  and  n 

are  dummy  variables 

V 

<&.);  1?; 

Jp  ^ 

From  the  definition  of  proportional  loading 

(111. 10) 

(111. 11) 


(III. 12) 


(III. 13) 


(III. 14) 
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Then  the  definition  of  the  current  stiffness  parameter  is 


e)  f  )  • 

~rr 

(pV 


(III. 15) 


Since  this  research  is  based  on  finite  element  analysis,  the  tensors 
a .  and  are  of  finite  dimension  n  ,  where  is  the  number  of  degrees 

of  freedom.  Since  these  tensors  are  finite  dimensional,  they  can  be 
represented  as  vectors.  In  vector  form,  the  current  stiffness  parameter  is 


dj~ 

<>ur 

jp 


(III. 16) 


The  stability  state  of  the  structure  can  be  evaluated  by  examining  the 
sign  of  the  current  stiffness  parameter  as  previously  explained.  This 
parameter  can  only  be  used  when  proportional  loading  is  used. 
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